;************************************************
; NCL Graphics: 
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"

;************************************************
begin

  load "load_setup.ncl"

  vname  = "INFrehet" 
  vunits = "" 
  vscale = 1.

  pa = basepa+ca+"/"
  pb = basepb+cb+"/"
  fa = ca+ma
  fb = cb+mb
  casa  = (/"A: "+na, "B: "+nb/)
  caseM  = (/pa+fa,pb+fb/)
  case  = (/caseM(0),caseM(1)/)
  Ncase = dimsizes(case)
  plotname = "figure_"+vname+"_zonal_"+shortname

  print(" ")
  print(" ")
  print("..................................................................")
  print(" ")
  print(" ")
  print(" varname   : "+plotname)
  print(" ctrl      : "+case(0))
  print(" test      : "+case(1))
  print(" plot name : "+plotname)
  print(" ")
  print(" ")

;************************************************
; Open netCDF file and read in dimensions
;************************************************

   filename=case(0)+".nc"
   fi=addfile(filename,"r")
   lev   = fi->lev
   lat   = fi->lat
   lon   = fi->lon 
   nlev  = dimsizes(lev)                           
   nlat  = dimsizes(lat)
   nlon  = dimsizes(lon)
   ;print("nlev "+sprinti("%7.2i",nlev))
   ;print("nlat "+sprinti("%7.2i",nlat))
   ;print("nlon "+sprinti("%7.2i",nlon))
   
;************************************************
; set default res
; res1 for gsn_csm_contour
; resp for gsn_panel
;************************************************
  wks   = gsn_open_wks (plotfmt,plotname)             ; open workstation
; gsn_define_colormap(wks,"gui_default")          ; choose colormap
  gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")          ; choose colormap
  res1                      = True                ; plot mods desired
  res1@gsnFrame             = False               ; advance page
  res1@gsnDraw              = False                
  res1@gsnSpreadColors      = True                 ; use full colormap 
  res1@cnFillOn             = True                 ; turn on color
  res1@cnLinesOn            = False                ; turn on/off contour lines
  res1@cnFillMode           = "AreaFill"           ; default "AreaFill" "RasterFill"
  res1@lbLabelBarOn          = True               ; turn on/off the labelbar
  res1@lbOrientation         = "Vertical"        ; labelbar Orientation "Vertical"
  res1@lbLabelsOn            = True               ; turn on/off the label of labelbar
  res1@lbBoxMinorExtentF     = 0.40                ; area allotted to each box of the LabelBar    
  res1@lbLabelFontHeightF    = 0.02                                               
  res1@lbLabelAutoStride     = False                ; optimal labels         
  res1@cnLineLabelsOn        = False                ; turn on line labels
  res1@cnLevelSelectionMode = "ExplicitLevels"      ; "ManualLevels""ExplicitLevels"
  ;;res1@cnLevels             = (/20,50,100,150,200,300,500,700,1000/)
  res1@cnLevels             = (/0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7/) 
  res1@cnLevelFlags         = False
  res1@tiXAxisFontHeightF    = 0.020               ; XAxis Font
  res1@tiYAxisFontHeightF    = 0.018               ; YAxis Font
  res1@tiYAxisString         = "Approximate Pressure(hPa)"
  ;res1@tiXAxisString         = "Latitude"
  res1@vpXF                  = 0.12     ; x location left
  res1@vpYF                  = 0.88     ; y location  up
  res1@vpWidthF              = 0.65     ; width
  res1@vpHeightF             = 0.35     ; height
  res1@gsnStringFontHeightF  = 0.02
  res1@tiMainString          = " " 
  res1@tiMainFontHeightF     = 0.03
  res1@tiMainFontThicknessF  = 0.35
  res1@trYReverse            = True
  res1@gsnYAxisIrregular2Linear = True

  resP                 = True          
  resP@gsnPanelYWhiteSpacePercent = 1  ; draw panel with white space added
  resP@gsnPanelXWhiteSpacePercent = 5  ; draw panel with white space added
  resP@txString            = vname 
  resP@gsnPanelLabelBar    = False                ; add common colorbar
  resP@lbLabelFontHeightF  = 0.012               ; labels size
  resP@lbBoxMinorExtentF   = 0.30                ; area allotted to each box of the LabelBar 
  
;************************************************
; read in variables & create plot
;************************************************

  plot = new(Ncase+2,graphic)


  res1@gsnSpreadColors      = True                 ; use full colormap 
  res1@gsnSpreadColorStart  = 20                    ; start at color 10
  res1@gsnSpreadColorEnd    = 195                    ; end at color 96
  
  do i=0,Ncase-1
     filename=case(i)+".nc"
     fi=addfile(filename,"r")
     ;print(" read "+case(i))
     VarLevLattema = dim_avg_Wrap(fi->$vname$(time|:,lev|:, lat|:, lon|:))
     VarLevLattemp = VarLevLattema(0,:,:) 
     VarLevLattemp = where(VarLevLattemp.gt.0.,VarLevLattemp,-999) 
     VarLevLattemp@_FillValue = -999
     VarLevLattemp=VarLevLattemp*vscale
     res1@gsnLeftString = casa(i)
     res1@gsnRightString= vunits                
     plot(i) = gsn_csm_contour(wks,VarLevLattemp({lev|0:1000},{lat|-90:90}), res1) 
     if(i.eq.0) then
        vara = VarLevLattemp
     end if
     if(i.eq.1) then
        varb = VarLevLattemp
     end if
     delete(VarLevLattemp)
  end do

  vara = where(vara.gt.0.,vara,-999) 
  vara@_FillValue = -999

  varc = vara
  varc = 100.*(varb - vara)/vara
  vard = vara
  vard = varb - vara

  ;;varc = where(vara.gt.5.,varc,-999) 

  res1@gsnLeftString = "(B-A)/A"
  res1@gsnRightString= "%"
  delete(res1@cnLevels)
  res1@cnLevels             = (/-40,-20,-10,-5,5,10,20,40/) ;;,60,100,200/)
  plot(2) = gsn_csm_contour(wks,varc({lev|0:1000},{lat|-90:90}), res1)

  res1@gsnLeftString = "B-A"
  res1@gsnRightString= vunits
  delete(res1@cnLevels)
  ;;res1@cnLevels             = (/-50,-20,-10,10,20,50,100,200,500,1000,2000/)
  res1@cnLevels             = 1.e-3 * (/-20,-10,-5,-2,-1,1,2,5,10,20/) 
  plot(3) = gsn_csm_contour(wks,vard({lev|0:1000},{lat|-90:90}), res1)

  
  gsn_panel(wks,plot,(/2,2/),resP)
  
  ;print("ps file : "+vname)
  
end


